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Abstract 



We present numerical calculations of the Local Density of Optical States (LDOS) in the near field of disor- 
dered plasmonic films. The calculations are based on an integral volume method, that takes into account 
polarization and retardation effects, and allows us to discriminate radiative and non-radiative contributions 
to the LDOS. At short distance, the LDOS fluctuations are dominated by non-radiative channels, showing 
that changes in the spontaneous dynamics of dipole emitters are driven by non-radiative coupling to plas- 
mon modes. Maps of radiative and non-radiative LDOS exhibit strong fluctuations, but with substantially 
different spatial distributions. 
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1. Introduction 

Disordered plasmonic films obtained by evapo- 
rating noble metals on a substrate are known to 
exhibit unusual optical properties Close to the 
percolation threshold, metallic clusters with fractal 
perimeters leads to the emergence of subwavelcngth 
areas supporting enhanced electric field, commonly 
called hot spots [Ij. These randomly distributed 
localized fields turned out to be very promising for 
sensing 0, Q] , subwavelength focusing Q , or non- 
linear optics [6[. Although several theoretical and 
numerical works have been reported on the subject, 
the question of the Local Density of Optical States 
(LDOS) has been hardly addressed. 

It has been known for long that the decay rate 
of a fluorescent emitter depends on its electromag- 
netic environment 0, 01 1 the dependence being de- 
scribed by the LDOS p(r,a;), with r the location 
of the emitter and uj the emission frequency. In- 
deed, the lifetime t of the excited state of a dipole 
emitter with transition dipole p is given in pertur- 
bation theory by 1/t = 7rw|pp/9(r, cj)/(3eofi.) where 
eo is the vacuum permittivity and h the reduced 



Planck constant. Thus the LDOS can be directly 
probed experimentally by measuring t. In a dis- 
ordered medium, changes in the LDOS probe the 
local environment 0, [lO, [Hi , the photon transport 
regime [H, [l^ or drive long-range correlations of 
speckle patterns [3, [Hi- Recently, LDOS statis- 
tics in the vicinity of disordered films have been 
studied experimentally [l^. Enhanced LDOS fluc- 
tuations have been observed close to the percolation 
threshold, in a regime where the film morphology 
is controlled by fractal clusters. These enhanced 
fluctuations have been qualitatively associated to 
localized plasmon modes. Theoretical and numer- 
ical studies of semi-continuous disordered metallic 
fllms are very often based on approximations, such 
as mean-field theories |17[ or quasi-static calcula- 
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tions [3] ■ An exact numerical approach has been 
reported recently using a FDTD (Finite-Difference 
Time-Domain) scheme [lij . 

In this paper, we present numerical calculations 
of the LDOS in the vicinity of disordered metallic 
films based on an integral volume method. This 
exact formulation is limited only by the discretiza- 
tion of the films into finite size cells. The numerical 
algorithm is divided into two steps. Firstly, we use 
a Monte-Carlo algorithm to simulate the growth of 
a gold film under an evaporation/deposition pro- 
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cess, and check that the geometrical properties of 
the film near the percolation threshold are in good 
agreement with experimental observations. Sec- 
ondly, we solve Maxwell's equations in 3D, taking 
into account polarization and retardation effects, 
which allows us to compute maps and statistical 
distributions of the LDOS. The computations are 
in agreement with known experimental results. The 
approach allows us to split the LDOS into its radia- 
tive and non-radiative contributions, and to discuss 
their relative contributions to the spatial fluctua- 
tions of the LDOS, which is the main focus of this 
work. 

2. Numerical approach 

2.1. Generation of disordered films 

Our flrst goal is to generate numerically disor- 
dered metallic films that have the same properties 
as the experimental evaporated metallic films. To 
do so, we use a kinetic Monte-Carlo algorithm, as 
proposed in The idea is to randomly deposit 
5-nm large gold particles on a square grid via an 
iterative algorithm, and let the particles diffuse un- 
der the influence of an interaction potential until 
a stable geometry is reached. At every iteration of 
the algorithm, we randomly choose either to deposit 
a new particle (probability pq) or to make a parti- 
cle on the grid jump to a more stable neighbour 
site (probability Pi^j to scatter from site i to site 
j). Using the normalization po + J2i,j^iPi~^j = 
we only need to pick a random number out of [0, 1] 
to determine the relative weight of each process. 
More precisely, the probability to deposit a par- 
ticle reads po = NF, where N is the number of 
particles that remains to be deposited in order to 
reach the prescribed filling fraction, and is a con- 
stant (with dimension s~^) modeling the experi- 
mental deposition rate. The probability for a par- 
ticle located on site i to jump to the neighbour 
site j reads p, 



eyip[—AEi-^j/{kBT)], where 
ks is the Boltzmann constant, T the temperature 
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the activation energy 
C omp uting AEi^j is a complex issue for 
22| . and is not possible from flrst prin- 
ciples for nanometer size particles. In the present 
approach, we have chosen to deal with a rescaled 
atomic potential that renormalizes the energy bar- 
rier in order to apply to a nanoparticle. We assume 
that AEi^j = a{Ei — Ej), where a is a positive 
dimensionless adjustable parameter taking into ac- 
count the influence of the substrate and the scaling. 



Ei is the rescaled "atomic" potential of a particle 
located on site i, which is allowed to jump to the 
neighbour site j if Ei > Ej . This potential is given 
by the following expression based on a tight-binding 
second moment method 12311: 



E, = A^exp[-p(ry/ro - 1)] 



1/2 



B^^exp[-29(r,,/ro-l)] ^ . (1) 

In this expression, rg is the size of one particle, 
rij the distance between two sites i and j and 
A, B, p and q are constants that were tabulated 
for atoms [23[. The iterative deposition process 



is stopped when all particles have been deposited 
(so that the prescribed fllling fraction has been 
reached) and no particle can move to a more stable 
site. 

Three examples of fllms, with a lateral size of 
375 nm and three different surface fllling fractions 
/, are shown in Fig. [TJ When the fllling fraction 



> ■< * * 
ft', _^ 



/ = 20 % 




375 nm 



f = 50% 



Figure 1: Numerically generated gold films for three different 
filling fractions / (gold is represented in dark). The param- 
eters for the computation are: T = 300 K, a = 2. 58.10"^, 
F = 10i*s-i, A = 0.2061 eV, B = 1.79 eV, p = 10.229, 
q = 4.036. 



increases, a continuous metallic path appears link- 
ing two sides of the sample (percolation). A very 
important feature of the disordered metallic fllms 
is the apparition of clusters with fractal perimeter 
near the percolation threshold [![. The perimeter 
P of a cluster is said to be fractal when Pftactai oc 
S^/^, where S is the cluster surface and D is a non- 
integer number called fractal dimension . Usual 
euelidian 2D surfaces have a perimeter satisfying 



Pc 



euelidian 



OC S^^^. It has been shown experimentally 
that on disordered metallic fllms, the fractal dimen- 
sion is D = 1.88 [2^. To check this feature, we 
generated 100 fllms with fllling fractions / = 20 % 
and / = 50 %. We extracted the perimeter and sur- 
face of all clusters in all numerically generated fllms. 
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Wc show in Fig. [2] the location of each cluster in a 
perimeter/surface diagram, in a log-log scale (each 
blue cross corresponds to one cluster) , for both fill- 
ing fractions. For low filling-fraction, every cluster 
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Figure 2: Distribution in a pcrimctcr/surfacc diagram of 
the clusters taken out from 100 numerically generated films. 
Left: filling fraction / = 20%. Right: filling fraction / = 
50 % . The red solid line and green dotted line are guides for 
the eye, corresponding to P = 75^/^ and P = 0.28S'^'**/2 ^ 
respectively. 



has an cuclidian perimeter {D = 1). For filling frac- 
tion / = 50 %, we clearly see the existence of fractal 
clusters with D ~ 1.88. This result, already shown 
in [2^, is a strong evidence that the geometrical 
features of experimental disordered films are well 
described by the numerical generation method. 

2.2. Expression of the LDOS 

In order to compute the electric field and the 
LDOS on disordered films, we consider that a unit 
pixel of a numerically generated film as that shown 
in Fig. [1] is a 5-nm size gold cube described by its 
dielectric constant e(w), taken from [2^. To com- 
pute the LDOS p(ro,a;), we have to compute the 
imaginary part of the dyadic electric Green func- 
tion G at the position of the emitter [23| • The 
normalized LDOS reads: 



P_ 

Po 



fco 



Im[TrG(ro,ro,cj)] 



(2) 



fco = w/c = 2ttc/X and po = lj /{-k c^) is the LDOS 
in free space. The dyadic Green function G con- 
nects an electric dipole p at position r' to the radi- 
ated electric field at position r through the relation 
E(r, oj) = /xoa;^G(r, r', a;)p. It describes the elec- 
tromagnetic response of the environment. 

2.3. Calculation of the dyadic Green function 

To compute the dyadic Green function G in the 
presence of the film, we consider an electric dipole 



p located at position rg , and use a volume integral 
method [28| . The electric field at any point r obeys 
the volume integral equation (Lippmann-Schwinger 
equation) 

E(r, uj) = fiouj^Go{r, Tq, uj)p 

+ kl[e{u)^l]( Go(r,r',c^)E(r',c^)dV (3) 
Jv 

where V is the volume occupied by the metallic 
film. Go is the dyadic Green function of free space, 
given by [ii [13 



Go(r,r',c^) = PV 



VV 

'''0 



exp(ifco-R) 



^<'- '''si- 



ll) 



where R= | r — r' | , P V denotes the Principal Value 
operator and 5 is the Dirac delta function. In or- 
der to solve the integral equation numerically, we 
discretize V into cells of size A, and assume that 
the electric field is constant in each cell (the vol- 
ume of cell number j will be denoted by Vj). For 
all calculations presented in this paper, A is set 
to 2.5 nm so that each gold cube is divided into 
eight cells. To improve convergence of the numer- 
ical computation, we integrate the Green dyadic 
on the cell volume (moment method) and define 
Gjf = L Go(rj,r',a;)d3r'. To calculate the elec- 
trie field in each cell, we have to solve the following 
linear system: 

{I - kl [e{u) 1] Gjf } E^-kl Hu;) - 1] ^ G'f E, 

= ^o'^^Go(i'i,ro,tj)p. (5) 

The solution leads to the expression of the three 
components of the electric field E.^ in cell number 
i, for all i. The computation of GJ"' has to be 
performed with care, due to the singularity of the 
Green function Go at the origin. This can be done 
in Fourie r sp ace, using the Weyl expansion as ex- 
posed in [31[. Solving Eqs. ^ for three orthogonal 
orientations of the source dipole gives direct access 
to the full dyadic Green function G(r,ro,a;). The 
LDOS is deduced from Eq. ([2]) (note that the imag- 
inary part of the Green dyadic is not singular at the 
origin for ro in vacuum). 

The numerical approach also allows us to calcu- 
late separately the radiative LDOS pR (which is 
proportional to the far-field power radiated by the 
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dipole source) and the non-radiative LDOS pnr 
(which is proportional to the power absorbed inside 
the metal) [s^ [s^ . Energy conservation requires 
that p — pr + pnr, so that only two quantities need 
to be calculated. We can compute the normalized 
non-radiative LDOS from 

^ = 11^ /.i^'-' 

and then deduce by subtraction. Equation ([5]) 
is discretized the same way as Eq. . ^'^From such 
calculations it is possible to address the contribu- 
tion of radiative and non-radiative modes to the 
LDOS, as we will see. This is an important issue 
in the understanding of the optical response of dis- 
ordered metallic films, and their use for the control 
of the dynamics of fluorescent sources. 

3. Results 



on subwavelength areas. The existence of local en- 
hancements of the electric field intensity (hot spots) 
is a well-known result, that was observed before in 
experiments 0] . These local field enhancements di- 
rectly translate into LDOS enhancements, leading 
to strongly fluctuating LDOS patterns. Another 
interesting output of the calculations is that at a 
distance 40 nm above a film with 375 nm lateral 
extension, LDOS spatial fluctuations are mainly 
due to non-radiative channels (this can be seen by 
comparing the standard deviations of pnr and pR 
in Fig. [3]). Moreover, the spatial distribution of 
the radiative LDOS pR is completely different from 
that of the non-radiative LDOS pnr- In a fluores- 
cence experiment using single nanoscale emitters, 
this means that the trade-off between radiative and 
non-radiative decay is dependent on the emitter po- 
sition. The apparent quantum yield also becomes a 
spatially fluctuating quantity, with expected strong 
fluctuations. 



3.1. Mapping the LDOS and its radiative and non- 
radiative components 

Using the approach described in section [21 we 
computed maps of the total, radiative and non- 
radiative LDOS at 40 nm above numerically gen- 
erated disordered metallic films. This distance has 
been chosen since it provides substantial near-field 
effects and remains compatible with standard com- 
putational resources. A full study of the distance 
dependence, both theoretical and experimental, will 
be published elsewhere 34|. The results are shown 
in Fig. [31 We clearly see that near the percola- 
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Figure 3; Maps of the total (p), non-radiative (pnr) a-iid 
radiative LDOS (pr) normalized by the LDOS in vaeuum 
(po) at 40 nm distanee above two films with filling fraetions 
/ = 20% and / = 50%. The wavelength is A = 780 nm. 
Note that the color scale is different for every map. 



tion threshold (film with / = 50%), complicated 
LDOS structures appear, with local enhancements 



3.2. Statistical distributions of p, pNR o,nd /9r 

The existence of localized modes on disordered 
metallic films was recently studied experimen- 
tally measuring the statistical distribution of the 
LDOS It was shown that the apparition of 

fractal clusters was correlated to enhanced fluctua- 
tions of the LDOS, that are a direct signature of the 
presence of spatially localized held distributions. 
We computed the statistical distribution of the to- 
tal, non-radiative and radiative LDOS, for two col- 
lections of films of lateral size 375 nm with filling 
fractions / = 20% and / 50%. For each filling 
fraction, we generated 60 different films and com- 
puted the value of the LDOS at a distance 40 nm 
above the center of the film. The histograms are 
shown in Fig. [31 ^From the calculations, we re- 
cover the enhanced fiuctuations of the LDOS ob- 
served in [l^l close to the percolation threshold. A 
comparison of the histograms for p, pnr and pR 
also confirms that at a distance 40 nm from the 
film, the LDOS fluctuations are mainly driven by 
non-radiative channels, as already discussed in sec- 
tion [3TT1 Finally we note that the computations are 
performed on samples with lateral size on the or- 
der of A/2, so that the LDOS spatial distribution 
might be affected by finite-size effects. Although 
not shown for brevity, we have performed compu- 
tations with sample sizes from f50nm to 375 nm. 
These computations have shown that although the 
statistical distribution of pR is size-dependent in 
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Figure 4: Histograms of the total (p), non-radiative (pnr) 
and radiative LDOS (pr) normalized by the LDOS in vac- 
uum (po ) at 40 nm distance above two series of films of same 
filling fraction (red: / = 20%; blue: / = 50%). Every 
generated film has a lateral size of 375 nm. 



this regime, the distribution of p and pnr are quite 
robust. 

3.3. Correlation between LDOS hot spots and film 
topography 

To get more insight about the origin of the lo- 
cahzed LDOS (or intensity) enhancements, we su- 
perimpose the maps of the total normahzed LDOS 



and the topography of the films, as shown in Fig.[SJ 
The maps clearly show that at low filling fraction 
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Figure 5: Maps of the total normalized LDOS at a distance 
40 nm represented on top of the film topography (gold is 
represented with black color). Wavelength A = 780 nm. Left: 
/ = 20%. Right: / = 50%. 



(left), classical plasmon resonances of isolated par- 
ticles are responsible for local enhancements of the 
LDOS. Near the percolation threshold (right), the 
origin of the LDOS structure is more complex. The 
non-trivial relation between the topography and the 
location of localized field enhancements is sustained 
by collective interactions. Finding a simple model 
to understand this connection is still an open issue. 

4. Conclusion 

In conclusion, we have presented exact 3D nu- 
merical calculations of maps and statistical distri- 
butions of the LDOS in the near-field of disor- 
dered plasmonic films. The calculations describe 
the well-known existence of localized enhancements 
of the near-field intensity and the LDOS on sub- 
wavelength areas, for filling fractions close to the 
percolation threshold. The method also permits a 
calculation of the radiative and non-radiative con- 
tributions to the LDOS. We have shown that at 
a distance 40 nm above the film (near- field zone), 
the LDOS fluctuations are chiefly driven by non- 
radiative channels. Nevertheless, both radiative 
and non-radiative LDOS exhibit strong spatial fluc- 
tuations, with completely different spatial distribu- 
tions. Understanding the trade-off between radia- 
tive and non-radiative channels is a key issue for 
the understanding of the optical properties of dis- 
ordered plasmonic films, and their use as sensors, 
absorbers or new materials for the control of light 
emission. 

References 

[f] V. M. Shalaev, Nonlinear Optics of Random Me- 
dia: Fractal Composites and Metal-Dielectric Films, 



5 



Springer Tracts in Modern Physics, Berlin Heidelberg 
(2000). 

[2] S. Gresillon, L. Aigouy, A. C. Boccara, J. C. Rivoal, X. 
Quelin, C. Desmarest, P. Gadenne, V. A. Shubin, A. 
K. Sarychev, and V. M. Shalaev, Phys. Rev. Lett. 82, 
4520 (1999). 

[3] M. Moskovits, Rev. Mod. Phys. 57, 783 (1985). 
[4] E. Fort and S. Gresillon, J. Phys. D: Appl. Phys. 41, 
013001 (2008). 

[5] X. Li and M.I. Stockman, Phys. Rev. B, 77 195109 
(2008). 

[6] A. K. Sarychev and V. M. Shalaev, Phys. Rep. 335, 
275 (2000). 

[7] E. M. Purcell, Phys. Rev. 69, 681 (1946). 

[8] R. R. Chance, A. Prock, and R. Sylbey, Adv. Chem. 

Phys. 37, 1 (1978). 
[9] L. S. Froufe-Perez, R. Carminati and J. J. Saenz, Phys. 

Rev. A 76, 013835 (2007). 
[10] L. S. Froufe-Perez and R. Carminati, Phys. Status So- 

lidi A 205, 1258 (2008). 
[11] R. Sapienza, P. Bondareff, R. Pierrat, B. Habert, R. 

Carminati and N. F. van Hulst, Phys. Rev. Lett. 106, 

163902 (2011). 

[12] H. Schomerus, M. Titov, P. W. Brouwer, C. W. J. 

Beenakker, Phys. Rev. B 65, 121101 (2002). 
[13] R. Pierrat and R. Carminati, Phys. Rev. A 81, 063802 

(2010). 

[14] B. A. van Tiggelen and S. E. Skipetrov, Phys. Rev. E 

73, 045601(R) (2006). 
[15] A. Cazc, R. Pierrat and R. Carminati, Phys. Rev. A 

82, 043823 (2010). 
[16] V. Krachmalnicoff, E. Castanio, Y. De Wilde and R. 

Carminati, Phys. Rev. Lett. 105, 183901 (2010). 
[17] X. C. Zeng, D. J. Bergman, P. M. Hui, D. Stroud, Phys. 

Rev. B 38, 10970 (1988). 
[18] M. I. Stockman, S. V. Faleev, and D. J. Bergman, Phys. 

Rev. Lett. 87, 167401 (2001). 
[19] U. K. Chettiar, P. Nyga, M. D. Thoreson, A. V. Kildi- 

shev, V. P. Drachev, V. M. Shalaev, Appl. Phys. B 100, 

159 (2010). 

[20] J. Aubineau, PhD thesis, University of Versailles Saint- 

Qucntin (2005). 
[21] R. Fcrrando and G. Trcglia, Phys. Rev. B 50, 12104 

(1994). 

[22] C. Mottet, R. Fcrrando, F. Hontinfinde and A.C. Levi, 
Surf. Sci. 417, 220 (1998). 

[23] F. Clcri and V. Rosato, Phys. Rev. B 48, 22-33 (1993). 

[24] B. Mandelbrot, The Fractal Geometry of Nature, Free- 
man, New York, USA (1983). 

[25] P. Gadenne, PhD thesis, University of Paris VI (1986). 

[26] E. D. Palik, Handbook of Optical constants of solids. 
Academic Press, San Diego, USA (1998). 

[27] J. M. Wylie and J. E. Sipe, Phys. Rev. A 30, 1185 
(1984). 

[28] R.F. Harrington, Field Computations by Moment 

Methods, IEEE Press, New York, USA (1992). 
[29] J. van Bladel, Singular Electromagnetic Fields and 

Sources, Clarendon, Oxford (1991). 
[30] A.D. Yaghjian, Proc. IEEE 68, 248(1980). 
[31] P. C. Chaumet, A. Scntonac and A. Rahmani, 

Phys. Rev. E 70, 036606 (2004). 
[32] R.X. Bian, R.C. Dunn, X. Sunney Xie, P.T. Leung, 

Phys. Rev. Lett. 75, 4772 (1995). 
[33] V.V. Klimov, M. Ducloy, V.S. Letokhov, Quant. Electr. 

31, 569 (2001). 



[34] E. Castanie, V. Krachmalnicoff, A. Cazc, R. Pierrat, Y. 
De Wilde and R. Carminati, submitted for publication 
( |http://arxiv.org/abs/1112.13T0| |. 



6 



